Thermocapillary migration and interactions of two nondeformable droplets 
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A numerical study on interactions ol two spherical drops in thermocapillary migration in mi- 
crogravity is presented. Finite-difference methods were adopted and the interfaces of drops were 
captured by the front-tracking technique. It is found that the arrangement of drops directly influ- 
ences their migrations and interaction, and that the motion of one drop is mainly determined by 
the disturbed temperature field because of the existence of the other drop. 
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Introduction 

Under the microgravity condition, the thermocapillary 
migration of droplets or bubbles in matrix liquid is caused 
by the nonuniform interface tension introduced by the 
temperature gradient. This motion is of great importance 
in material processing and other applications in space. 
The original work in this field was performed by Young et 
al. [l| ■ In their study, the inertial convection and thermal 
convection are neglected (the so-called YGB Model), and 
the derived migration velocity is 



Vygb — 



2U 



(2 + 3fx d /n b )(2 + k d /k b ) 



(1) 



Here, U is the reference velocity defined by the bal- 
ance of thermocapillary force and viscosity force on the 
drop/bubble: 

U = |or || VT^a//^, 

fi is the kinematic viscosity, k the thermal conductivity, 
<tt the rate of change of interfacial tension with temper- 
ature, VTqo the temperature gradient imposed on ma- 
trix liquid, and a the radius of the drop or bubble. The 
symbols with the subscript d mean the parameters of the 
droplet /bubble, and those of the bulk liquid are indicated 
by the subscript b. 

After YGB, there are many other studies on the ther- 
mocapillary motion of isolated drop/bubble (see [2| and 
references therein). In practice, it is common to have 
two or more drops/bubbles in the continues phase, so it 
is necessary to study their interactions. The first axisym- 
metric investigation of two thermocapillary bubbles was 



conducted by Meyyappan et a/.0], using the bipolar co- 
ordinate. It was found that the smaller bubble always 
moves faster than the isolated drop while the bigger one 
moves slightly slower. Meyyappan and SubramanianQ 
extended the above work to arbitrarily placed bubbles. 
Balasubramaniam and SubramanianQ assumed that two 
bubbles migrated in the potential flow (namely, the re- 
lated Re number is very large), and the matched asymp- 
totic analysis was adopted to solve the energy equation 
with large Ma numbers. It was found that the thermal 
wake of the leading bubble will disturb the temperature 
field around the trailing bubble and reduce its velocity. 
Interactions between two spherical droplets were firstly 
studied by Anderson with a reflection method^. It was 
found that interactions between droplets driven by ther- 
mocapillary effects are much weaker than those of sed- 
imentation. Ken and ChenQ analyzed the axisymmet- 
ric motion of two droplets in the bishperical coordinate, 
and their later combined analytical-numerical study was 
about a finite chain of spherical droplets along the line of 
their centers Q . Interactions of two deformable droplets 
in the axisymmetric coordinate were studied by Zhou 
and Davis [9|. Thermocapillary interactions of droplets 
or bubbles toward a hot wall at finite Reynolds and 
Marangoni numbers were numerically studied by Nas et 
al. 

0,1 It was found that bubbles and light drops 
line up perpendicular to the temperature gradient and 
are evenly spaced in the horizontal direction. A space 
experiment observed that a small leading drop could re- 
tard the movement of the big trailing drop 
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So far as we know, there is no systematic study on the 
interaction of two arbitrarily placed drops in the thermo- 
capillary research, and it will be the main subject here. 
We focus our investigation on two droplets with equal 
size and the same physical parameters (kinematic viscos- 
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ity, thermal diffusivity, density, and specific heat). The 
governing equations and numerical methods will be intro- 
duced in the next section, detailed numerical models and 
parameters of simulations are in section 2, and the results 
when the inertia and thermal convection are ignorablc or 
not, are discussed in section 3 and 4, respectively. 



I. GOVERNING EQUATIONS AND 
NUMERICAL METHODS 

In the thermocapillary motion, the two droplets with 
the same radius a are surrounded by the bulk fluid in a 
rectangular box ft = [xq,xi] x [yo,yi] x [20,21] (Fig- ID- 
The box is closed by no-slip walls. The direction of the 
temperature gradient is along the z axis, and x = is 
treated in the drop centers. The governing equations for 
this problem are: 



V • u = 0, 
d(pu) 



dt 
PC P ( 



dT 
~dt 



V • (puu) = -Vp + V ■ (m(Vu + V T u)) 
fu-VT) = V-(JfcVT). 



Here, u = (u,v,w), x = (x,y,z) G tt. F a is the body 
force term calculated by integrating the surface tension 
across the interface [2(. Except the different material pa- 
rameters for the drop phase and the bulk phase, the con- 
servative equations above are valid for both phases. We 
define the nondimensional quantities as: 



u = u/U, 
p = p/( Pb U 2 ), 
A = 



x = x/a, t = t/{—), 

T = T/(\VT 00 \a), p = p/ Pb , (2) 



(J, /fib, k = k/k b , Cp = C p /Cpb 1 
F a = F a a/(p 1 U 2 ), Re=Ua/v b , Ma = Ua/n b . 

Here, v b = p b /p b is the kinematic viscosity, and K b = 
k b /(p b C pb ) the thermal diffusivity of the matrix liquid. 
The nondimensional equations can be written as: 



V • u = 
d(pu) 



dt 



V • (puu) 



-Vp- 



-V.(MVu- 



V T u)) 



pC p (§ + u.VT) = ±-V.(Wf). 
The boundary conditions for velocities are: 
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it 



X — Xq .. 



X\ ^\x — Xq^X\ ^\x — Xq^X\ 0i 



y=vo:Vi ~ v \y=vo,vi 



\Z — Zq,Zi 



W 



= 0, 
0. 



(3) 

(4) 
(5) 

(6) 




FIG. 1: The Sketch of two drops in the thermocapillary mi- 
gration. $ is the angle between the temperature gradient and 
the center line of two drops. 

For energy equation, the Dirichlet boundary condition is 
adopted: 



T\x=x — Tq 



\y=vo 



The initial conditions are: 



T\x=x 1 — To 
T \y=yi = To ' 



■z, 



(7) 



T\ s =zi =T + zi, 



u\t=o 



T 



t=0 



V\t=o = W\t =0 
:f + Z. 



0. 



(8) 



In the following, symbols without overbars will be 
adopted to denote non-dimensional values. 

II. NUMERICAL MODELS AND PARAMETERS 

Fig. [T] indicates initial positions of two drops. The 
horizontal and vertical distances between two drops are 
dy and dz, respectively. The traditional definition of the 
non-dimensional distances arc defined as S y = dy/2 and 
S z = dz/2 [2j, and S y o and S z o denote the initial val- 
ues of S y and S z . 9 indicates the point on drop inter- 
face in the x — plane: 9 = is the front stagnation 
and 9 — ir or 9 = —ir the rear stagnation. Points in 
the clockwise direction from front stagnation are denoted 
with 9 > 0, otherwise, 9 < 0. In the full three dimen- 
sional model, the computation zone is set to be 6 x 9 x 24 
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FIG. 2: The steady-state thermocapillary migration velocities 
of two drops with different initial distances. The horizontal 
ordinate is S z o for $ = 0, and S y o for $ = tt/2. Here, Re = 
Ma = 10" 3 . 




y/(Ua) = I 




(b) 



FIG. 3: Streamlines in a reference frame attached to the 
drop/sphere, (a) Potential flow, (b) Stokes flow. 

III. THERMOCAPILLARY MIGRATION OF 
TWO DROPS WITH NEGLECTED INERTIA 
AND THERMAL CONVECTION 



on a grid of 60 x 90 x 240. The time steps are 10~ 6 for 
Re = Ma = 10~ 3 , and 10~ 3 for other RekMa values. 
To save the computing time, the axisymmetric model is 



adopted in the cases of $ = 0|22j. In axisymmetric sim- 
ulations, dropl in hotter region will be called the leading 
drop and drop2 in colder region the trailing drop, and S 
and So are adopted to replace S z and S z o in the three- 
dimensional model. The computing domain is 6 x 24 with 
the resolution of 128 x 512. The time steps are 5 x 10~ 7 
for the simulations of Re = Ma = 10~ 3 and 2 x 10~ 4 
for all other RekMa values. In this paper, all material 
parameters of two drops are assumed to be the same. 

Generally speaking, the non-dimensional thermocapil- 
lary migration velocities of drops are quite small (about 
0.1), and the usually defined non-dimensioan velocity 
(namely u = u/U ma xi where U max is the maximum ve- 
locity in the flow field instead of U = lorllVToola/ Mb) i s 
even lower. The influence of the Re number in the cur- 
rent study is trivial and can be inferred from the results 
of the isolated drop. The role of the Re number will not 
be discussed here (simply set Re = 1 if not specified), 
and we will concentrate on the influences of thermal con- 
vection and initial distance. 

To have a clear idea of the interaction of two thermo- 
capillary drops, it is necessary to compare it with that 
of the isolated drop. In the following, the velocity of the 
isolated drop will be denoted as Wi S0 . 



In this section, small Re and Ma numbers (10~ 3 ) are 
adopted in simulations to compare with some previous 
analytical results. Fig. [2] shows the final migration ve- 
locities for $ = & 7r/2. In the case of $ = 0, velocities 
of the two drops are very close, but faster than Wi SO (see 
the solid line in Fig. [2]). This phenomenon is also found 
in the previous analytical and numerical studies 0, 
In the case of $ = tt/2, our three dimensional simula- 
tions show that the drop velocities in y or x direction are 
neglectable. The z direction velocities (W) of both drops 
are still the same, but they are slower than Wi SO (the 
dash line in Fig. [2]). In both cases, the increase of the 
initial drop distance makes the final drop velocity closer 
to W iso . 

Similarly, there are lots of researches on the inter- 
actions between two rigid spheres 15- 1^. Using pre- 
vious analytical linear results, we compare the differ- 
ence between drops and rigid spheres. Assume the rigid 
sphere/drop has a constant velocity U, and the original 
point is in the center of the rigid sphere/drop. In the 
spherical coordinate (r,8,<f>), the potential flow, which 
describes the thermocapillary motion of the drop, is writ- 
ten asQ: 



v r 



U —j- cos ( 



V0 



l r V . 

--U -it sm( 



(9) 



Here, 8 = is the motion direction. The motion of the 
rigid sphere is described by the Stokes flow|19j: 
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FIG. 4: The migration velocities of the leading drop, the trailing drop, and the isolated drop with Re = 1 and So = 1.5. 



v r (r,6) = lu(3- - %)cos0 
2 r r A 



IV. THERMOCAPILLARY MIGRATION OF 
TWO DROPS WITH FINITE INERTIA AND 
THERMAL CONVECTION 

La Influence of thermal convection for the cases of 

$ = 



v 6 {r,9) = —U(Z- + -r) 



sin t 



(10) 



It is clear that, with the increase of r, the velocities 
for the rigid sphere and drop are decaying in the mag- 
nitude of 0(l/r 3 ) and 0(l/r), respectively. The veloc- 
ity perturbation in potential flow spreads in a relatively 
small region (Fig. 02[a)), and has obvious directionality 
because the liquid in the front of the moving drop will 
be supplied to the back of the drop. For Stokes flow, the 
velocity perturbation spreads in a fairly large region, and 
the surrounding liquid tries to move with the sphere(Fig. 

Mb))- 

The drag coefficients of both rigid spheres are always 
lower than that of the isolated sphere 2CJ, l2l| , which 



means two spheres will move faster than an isolated one. 
When the nonlinear effect is strong, the main interest in 
the interaction between rigid spheres is drag coefficients 
changed by the wake flow behind the leading body. How- 
ever, in the thermocapillary study, the disturbed temper- 
ature field is the most important, and we will concentrate 
on it in the following section. 



Firstly, we study the influence of the Ma number on 
the interaction between two droplets in thermocapillary 
motion with Sq = 1.5 and Re = 1. 

In the case of Ma = 1, the leading drop is faster than 
the trailing drop, but both of them move faster than the 
isolated drop (Fig. HJa)). This is similar to what we 
have discussed in the last section. For Ma = 20, the 
leading drop is faster than the isolated drop throughout 
the simulation, while the speed of the trailing drop is up 
to 8% lower than W lso (Fig. gfb)). For Ma = 100, the 
trailing drop is even slower, and its velocity is up to 20% 
lower than W iso (Fig. H^c)). 

Figs. [5] are the isotherms around the drops. When the 
Ma number is small, the isotherms around the drops are 
almost straight and evenly spaced throughout the simula- 
tion (Figs. EIa)(d)). When the Ma number is increased, 
isotherms near r — arch to the hotter region, and there 
is a closed cold zone arising in the droplet (Figs. 0Hc)(f)). 
The temperature at the rear stagnation point of the lead- 
ing drop is lower than that of the isolated drop, and the 
temperature gradient of the trailing drop is also reduced. 

Fig. |5] shows the temperature difference between the 
point on drop surface and the front stagnation (to get 
a clear idea of this difference, the corresponding value 
of the isolated drop is subtracted in Figs. El l8l and IM)) . 
It can be seen that the temperature difference between 
the front and rear stagnation points of the leading drop 
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(solid line) is larger than that of the isolated drop; on the 
contrary, the difference for the trailing drop (dashed line) 
is smaller than that of the isolated drop, and will decrease 
with increasing Ma numbers. Fig. [6] clearly shows that 
the influence of the thermal convection on trailing drops 
is much stronger than that on leading drops. 

When the heat convection is stronger, the influence on 
the trailing drop is also bigger. Hence, the separated 
distances (S — So) increase more rapidly for larger Ma 
numbers (Fig. [7|). Note that the separating speed of 
Ma = 100 before t = 60 is lower than that of Ma = 20. 
This is because the fluctuation process of migrating speed 
in the beginning is longer for larger Ma number [2|. For 
example, the thermal wake left by leading drop with 
Ma = 100 is not fully developed until t = 60 (Fig. 
[5fc)(f)), and the velocity difference between leading and 
trailing drops is not so large. 

Temperature differences for Ma — 20 at various mo- 
ments are shown in Fig. [8j At t = 20, the tempera- 
ture difference of the leading drop becomes bigger than 
that of isolated drop, and starts to become smaller af- 
terwards. When the thermal convection effect and the 
temperature disturbance caused by the leading drop are 
fully developed, the temperature difference of the trail- 
ing drop reaches its minimum at t = 40, and starts to 
become larger afterwards. 

B. Influence of initial distance for the cases of $ = 

Obviously, there will be no interaction between drops if 
they are far enough from each other, so it is interesting to 
know the critical initial distance at which both droplets 
migrate like an isolated one. 

Several different initial distances are tested for the case 
of Re = 1 and Ma = 20 (Figs. [9j. The curve of the 
leading drop with Sq — 3 is almost identical to that of 
the isolated drop (Fig. EJa)), while the critical initial 
distance is 5 for the trailing drop (Fig. EJb)). It seems 
that the thermal wake left by the leading drop affects 
the trailing drop at a longer distance than the distance 
at which the trailing drop interferes the leading drop. 

When the Ma number is increased to 100, the critical 
initial distance for the leading drop remains to be Sq = 3 
(Fig. [TOta)). On the other hand, the thermal wake left 
by the leading drop is much longer and the critical initial 
distance for the trailing drop seems longer than 5. 

Figs. [TTJ shows the time evolutions of distances be- 
tween two drops. It's clear that the smaller the initial 



distance between the two drops, the faster they will sep- 
arate from each other. In the later stage of simulations, 
it seems that the distance differences caused by various 
initial distances vanish, and the final distance between 
drops (of course, if they are not too far away apart) is 
determined by other parameters. 

C. Influence of thermal convection for the cases of 

$ / 

In this subsection, the full three dimensional problem 
with Re = 1, Ma — 20 and S y0 = S z0 = 1.35 is studied. 
In this simulation, both drops have zero velocities in the 
x direction. In the vertical direction, the upper dropl 
migrates slower than the isolated drop while the lower 
drop2 moves faster than the isolated drop (Fig. IT^T a)). 
As a result, the vertical distance between the two drops 
is always decreasing (see the solid line in Fig. [IS"]) . 

The isotherms at t = 60 are shown in Fig. [T2] It 
can be seen that the thermal convection of the lower left 
drop2 causes the bending of isotherms around the upper 
right dropl. Compared with the isolated drop, the tem- 
perature gradient on the lower part of the left drop2 is 
reduced, while that on the upper part of the right dropl 
is enlarged. In order to get a better understanding of the 
velocities in z direction (W), we study the temperature 
distributions at t = 60 in the x = plane (Figs. [Ml) . 
Compared with the isolated drop, the temperature dif- 
ference between front and rear stagnation points of drop2 
is larger, while that of dropl is smaller. 

The drop velocities in y direction (V) are plotted in 
Fig. Q2] (b). It is found that dropl is always trying to 
move away from drop2 in y direction. Drop2 moves to- 
wards dropl in the beginning, but starts to move away 
since t « 30. Because Vd r0 p2 is always larger than Vdropi, 
the horizontal distance between two drops is increasing 
throughout the simulation. 

D. Influence of the initial distance for the cases of 

$ / 

In this subsection, three sets of simulations with Re = 
1 and Ma — 20 starting from different initial distances 
are studied: 

1- (S y0 ,S z0 ) = (1.35,1.35); 

2. (S y0 ,S z0 ) = (1.1,1.35); 



r 

(d) 



r 

(e) 



r 

(f) 



FIG. 5: The isotherms around the two droplets with Re = 1 and So = 1.5. The first row is at t = 20, and the second row 
t = 60. The first column: Ma = 1; the second column: Ma = 20; the third column: Ma = 100. 



3. (S v o,S z0 ) = (1.35,1.25). 

It is clear that the smaller the initial horizontal distance, 
the larger the velocities of two drops in the y direction, 
and the bigger the temperature differences between the 
left and right sides of drops (T(0) - T{-6), Fig. Mb))- 

With different Re and Ma numbers, the evolutions of 
distances between two drops for easel are shown in Fig. 
[T5l It can be seen that the two droplets separate very 
slowly when Re = 1 and Ma = 1. When the Re num- 
ber is increased, the two drops get close faster in the 
vertical direction, while there are only trivial changes in 
separated distances for increasing Ma numbers. 

Generally speaking, in the z direction, the lower drop2 
moves faster than the isolated drop, while the upper 
dropl moves slower than the isolated drop, and thus S z is 
decreasing throughout any simulation in this subsection. 
If the simulation domain is big enough, drop2 would ex- 
ceed dropl in the z direction and slow down to a velocity 



smaller than that of dropl. Then dropl will start to catch 
up with drop2, and so on. Eventually, both droplets will 
reach a steady migration state when they are aligned hor- 
izontally, as indicated by Nas et al. 
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V. CONCLUSIONS 

In this paper, the interactions of two nondeformable 
droplets in thermocapillary motion are studied. When 
the inertia and thermal convections are neglected, two 
vertically-placed drops will move faster than the iso- 
lated drop, while two horizontally-placed drops will move 
slower. For the finite Ma number and <f> = 0, the lead- 
ing drop moves faster than the isolated drop, while the 
trailing drop migrates slower than the isolated drop due 
to the disturbed temperature field left by the leading 
drop. When the two drops are closer, their interaction is 
stronger, but this intensive interaction will not last long 
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FIG. 6: Temperature difference between the point 9 on in- 
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FIG. 7: Time evolutions of separation distances (S—So) 
between two drops for Re — 1 and So = 1.5. 
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FIG. 8: The time evolution of the temperature differ- 
ence for Re = 1, Ma = 20 and So = 1.5. 

because the velocity difference of two drops is also big. 
Once there is enough big gap between the two drops, 
they will migrate like the isolated drops. For the finite 
Ma number and <£> ^ 0, the motions of two droplets are 
still limited in the y — z plane. The upper dropl mi- 
grates slower while the lower drop2 migrates faster than 
the isolated drop, which results in a smaller vertical dis- 
tance and a bigger horizontal distance between the two 
drops. Tab. Q]sums up the velocities of dropl and drop2 
in the z direction {W\, W2) studied in this paper. 

Here, we only explore a few interacting mechanisms 
of two droplets with a limited number of parameters. 
A wider range of parameters as well as deeper physical 
explanations should be included in the future works. 
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FIG. 10: Re — 1, Ma = 100, time evolutions of drop velocities with various initial distances, (a) Leading drop, (b)trailing 
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FIG. 11: Time evolutions of the separating distances between drops, (a) Re — 1, Ma = 20, (b) Re = 1, Ma — 100. 
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FIG. 12: Re = 1, Ma — 20 and S y o = S z o = 1.35, evolution of migration velocities, (a) velocity in z direction and (b) velocity 
in y direction. 



16 




2 4 y 6 



FIG. 13: Isotherms at t = 60. Re = 1, Ma = 20 and S y0 = 
Szo = 1.35. 
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It has been realized that some nonaxisymmetrical behav- 
iors might arise for the isolated drop when the full three- 
dimensional simulating domain is fairly small [3]. In this pa- 
per, we adopt a fairly large domain, and assume the nonax- 
isymmetrical effect can be neglected for moderate parameters. 
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FIG. 14: The temperature distributing on the surface of drops in the x = plane at t = 60. Re = 1, Ma = 20, S y o = S z o = 1.35 
and &T{6)=T{6)-T iso {e). 




FIG. 15: The time evolution of the vertical and horizontal 
distances between two drops with S y o = S z o = 1.35. 
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FIG. 16: (a) Time evolutions of the vertical velocities of the drops with various initial distances, (b) the temperature differences 
between the left and right sides of the drops. 
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Parameters 


Re = Ma= 10~ 3 


Mo=l,fl£ = l 


Ma = 20, RE = 1 


Ma = 100, i?£ = 1 


rigid Spheres 


$ = 


Wi + ,W 2 + 


+ , W 2 + 


Wi+,W 2 - 


Wi+,W 2 - 


Wi+, W 2 + 


$ = 0.68, tt/4, 0.82 






W 1 -,W 2 + 




Wi+,1V 2 + 


$ = tt/2 


Wi-,W 2 - 








Wi+, w 2 + 



TABLE I: The velocities of drop 1 (Wi) and drop2 (W 2 ) are compared with Wiso. '+'/'—' means the velocity is bigger/smaller 
than Wiso- The velocities of rigid spheres in Stokes flow are listed in the last row, the bigger/smaller velocity stands for the 
smaller/bigger resistance than that on the isolated rigid sphere. 



